﻿using System;
using System.Diagnostics;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Nature.Numerics;
using System.Xml.Linq;

namespace Nature.Toolkit.GMix
{
    [DebuggerDisplay("{ComponentName} Nsp: {Nsp}, Nfr: {Nfr}")]
    public class CompositionObject : Nature.IMixtureDefinition
    {
        #region Private Fields
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        double[] m_spX;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        double[] m_spY;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        double m_molarMass;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        double[] m_frX;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        double[] m_frY;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        CompositionObject[] m_fractions;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        CHOFormula m_choFormula;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        private string m_speciesID;

        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        private Temperature? m_temperature;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        private Pressure? m_pressure;

        #endregion

        internal CompositionObject(ExpressionInfo expression, string speciesID)
        {
            this.Expression = expression;
            m_speciesID = speciesID;
            SetFractionsX(new double[0]);
            m_fractions = new CompositionObject[0];            
        }

        internal CompositionObject(ComponentInfo composition)
        {
            this.Component = composition;
            this.Expression = composition.Expression;            
        }

        internal MixProjection Projection { get; set; }

        internal void SetFractions(CompositionObject[] fractions)
        {
            m_fractions = fractions;
        }


        public string ComponentName
        {
            get 
            {
                return (Component != null 
                    ? Component.DisplayName
                    : string.Format("{0} - pure chemical species", m_speciesID));
            }
        }

        /// <summary>
        /// Number of fractions defining this composition
        /// </summary>
        public int Nfr
        {
            get { return m_fractions.Length; }
        }

        internal ComponentInfo Component { get; private set; }

        internal ExpressionInfo Expression { get; private set; }

        public double GetParameterValue(string parameterName)
        {
            if (this.Projection == null)
            {
                Debug.Assert(this.Component.IsRoot == false);
                throw new InvalidOperationException(
                        string.Format("Parameter values can be queried only on root components")
                    );
            }
            try
            {
                int parameterIndex = this.Expression.ParameterIndex(parameterName);
                return this.Projection.ParameterValues[parameterIndex];
            }
            catch (KeyNotFoundException)
            {
                var error = new StringBuilder();
                error.AppendFormat(@"""{0}"" is not recognized as a valid parameter name", parameterName ?? "null");
                if (this.Expression.ParameterNames.Any()) 
                {
                    error.AppendFormat(" The following parameters are found in the mixture expression:");
                    for (int i = 0; i < this.Expression.ParameterNames.Length; ++i)
                    {
                        error.AppendFormat("{0} {1}", i > 0 ? "," : string.Empty, this.Expression.ParameterNames[i]);
                    }
                }                
                throw new InvalidOperationException(error.ToString());
            }
        }

        /// <summary>
        /// Mole ractions array
        /// </summary>
        public ReadOnlyArray<double> X 
        {
            get { return m_spX; }
            internal set
            {
                var masses = Expression.DataSet.SpeciesMasses;
                m_spY = new double[masses.Length];
                double sum = 0.0;
                for (int isp = 0; isp < masses.Length; ++isp)
                {
                    double add = value[isp] * masses[isp];
                    sum += add;
                    m_spY[isp] = add;
                }
                for (int isp = 0; isp < masses.Length; ++isp) { m_spY[isp] /= sum; }
                m_spX = value.ToArray();
                m_molarMass = sum;
                OnCompositionReady();
            }
        }

        /// <summary>
        /// Mass ractions array
        /// </summary>
        public ReadOnlyArray<double> Y
        {
            get { return m_spY; }
            internal set
            {
                var masses = Expression.DataSet.SpeciesMasses;
                m_spX = new double[masses.Length];
                double sum = 0.0;
                for (int isp = 0; isp < masses.Length; ++isp)
                {
                    double add = value[isp] / masses[isp];
                    sum += add;
                    m_spX[isp] = add;
                }
                for (int isp = 0; isp < masses.Length; ++isp) { m_spX[isp] /= sum; }
                m_spY = value.ToArray();
                m_molarMass = 1.0 / sum;
                OnCompositionReady();
            }
        }

        public ReadOnlyArray<string> SpeciesIDs
        {
            get { return this.Expression.SpeciesIDs; }
        }

        public Temperature Temperature 
        {
            get { return m_temperature.Value; }
            internal set { m_temperature = value; }
        }

        public Pressure Pressure 
        {
            get { return m_pressure.Value; }
            internal set { m_pressure = value; }
        }

        public Temperature? BaseTemperature { get; internal set; }

        public Pressure? BasePressure { get; internal set; }

        /// <summary>
        /// Mole fractions of the components of the mixture
        /// </summary>
        public ReadOnlyArray<double> XComp { get { return m_frX; } }

        public ReadOnlyArray<double> YComp { get { return m_frY; } }

        public ReadOnlyArray<CompositionObject> Fractions { get { return m_fractions; } }

        public double MolarMass { get { return m_molarMass; } }

        public CHOFormula CHOFormula { get { return m_choFormula; } }

        public EquivalenceRatio? EquivalenceRatio { get; internal set; }

        public MixtureFraction? MixtureFraction { get; internal set; }

        public CompressionRatio? CompressionRatio { get; internal set; }

        public MachNumber? MachNumber { get; internal set; }

        /// <summary>
        /// Number of species
        /// </summary>
        public int Nsp 
        {
            get { return X.Length; }
        }


        internal void AddToY(double[] y, double fraction)
        {
            for (int isp = 0; isp < Expression.DataSet.Nsp; ++isp)
            {
                y[isp] += fraction * m_spY[isp];
            }
        }

        internal void AddToX(double[] x, double fraction)
        {
            for (int isp = 0; isp < Expression.DataSet.Nsp; ++isp)
            {
                x[isp] += fraction * m_spX[isp];
            }
        }

        internal void SetFractionsX(double[] x)
        {
            m_frX = new double[x.Length];
            m_frY = new double[x.Length];
            double mass = 0.0;
            for (int ifr = 0; ifr < x.Length; ++ifr)
            {
                double add = x[ifr] * m_fractions[ifr].MolarMass;
                m_frX[ifr] = x[ifr];
                m_frY[ifr] = add;
                mass += add;
            }
            for (int ifr = 0; ifr < x.Length; ++ifr) { m_frY[ifr] /= mass; }
        }

        internal void SetFractionsY(double[] y)
        {
            m_frX = new double[y.Length];
            m_frY = new double[y.Length];
            double mass = 0.0;
            for (int ifr = 0; ifr < y.Length; ++ifr)
            {
                double add = y[ifr] / m_fractions[ifr].MolarMass;
                m_frY[ifr] = y[ifr];
                m_frX[ifr] = add;
                mass += add;
            }
            for (int ifr = 0; ifr < y.Length; ++ifr) { m_frX[ifr] /= mass; }
        }
        


        private void OnCompositionReady()
        {
            double c = 0.0, h = 0.0, o = 0.0;
            for (int isp = 0; isp < Expression.DataSet.Nsp; ++isp )
            {
                var cho = Expression.SpeciesCHOFormulas[isp];
                double x = m_spX[isp];
                c += x * cho.Carbon;
                h += x * cho.Hydrogen;
                o += x * cho.Oxygen;
            }
            m_choFormula = new CHOFormula(c, h, o);
        }

        internal double CalcEntropy(double temperature)
        {
            using (var managedSpeciesArray = Expression.SpeciesArrayObjectPool.Acquire())
            {
                double[] s = managedSpeciesArray.Unwrap();
                Expression.CalcEntropies(s, temperature);
                double entropy = 0.0;
                for (int isp = 0; isp < Expression.Nsp; ++isp)
                {
                    double x = m_spX[isp];
                    entropy += x * (s[isp] - Math.Log(x + 1.0e-16));
                }
                return entropy * Constants.Rgas;
            }
        }


        internal double CalcEntropy(double temperature, double pressure)
        {
            double s = CalcEntropy(temperature);
            s -= Constants.Rgas * Math.Log(pressure / Constants.Atmosphere);
            return s;
        }

        internal double CalcEnthalpy(double temperature)
        {
            using (var arrayObject = Expression.SpeciesArrayObjectPool.Acquire())
            {
                double[] h = arrayObject.Unwrap();
                Expression.CalcEnthalpies(h, temperature);
                double enthalpy = 0.0;
                for (int isp = 0; isp < Expression.Nsp; ++isp)
                {
                    enthalpy += m_spX[isp] * h[isp];
                }
                return enthalpy * Constants.Rgas * temperature;
            }
        }


        internal Temperature CalcIsentropicCompressionTemperature(Pressure pressure)
        {
            if (this.Pressure == pressure)
            {
                return this.Temperature;
            }
            //using (var managedEntropySolver = Expression.ObjectPoolContainer.Acquire<EntropyConservationLawTemperatureSolver>())
            //{
            //    var solver = managedEntropySolver.Unwrap();
            //    solver.Function.Configure(this, this.Temperature, this.Pressure, pressure);
            //    double temperature = solver.Solve(this.Temperature);
            //    return new Temperature(temperature, Component.TemperatureUnits);
            //}

            using (var solverObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionSoneSolver>())
            using (var vectorObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionVector>())
            {
                var solver = solverObject.Unwrap();
                var vector = vectorObject.Unwrap();

                Temperature? baseTemperatureMemo = BaseTemperature;
                Pressure? basePressureMemo = BasePressure;
                try
                {
                    BaseTemperature = Temperature;
                    BasePressure = Pressure;
                    var ratios = solver.Function.Configure(this, pressure);
                    vector.CompressionRatio = ratios.DensityRatio;
                    vector.Pressure = pressure;
                    vector.Temperature = Temperature * ratios.TemperatureRatio;
                    object anchor = null;
                    solver.Solve(vector, ref anchor);
                    return new Temperature(vector.Temperature).Convert(Component.TemperatureUnits);
                }
                finally 
                {
                    BaseTemperature = baseTemperatureMemo;
                    BasePressure = basePressureMemo;
                }
            }
        }

        internal double CalcEnthalpy(Pressure pressure)
        {
            Temperature t = CalcIsentropicCompressionTemperature(pressure);
            return CalcEnthalpy(t);
        }


        internal void ApplyAdiabaticMixing(Pressure pressure)
        {
            using (var managedEnthalpySolver = Expression.ObjectPoolContainer.Acquire<EnthalpyConservationLawSolver>())
            {
                var enthalpySolver = managedEnthalpySolver.Unwrap();
                double fixedEnthalpy = 0.0;
                double tempEstimate = 0.0;
                for (int ifr = 0; ifr < Nfr; ++ifr)
                {
                    var fraction = m_fractions[ifr];
                    Temperature t = fraction.CalcIsentropicCompressionTemperature(pressure);
                    fixedEnthalpy += m_frX[ifr] * fraction.CalcEnthalpy(t);
                    tempEstimate += m_frX[ifr] * t;
                }

                enthalpySolver.Function.Configure(this, fixedEnthalpy);
                double temperature = enthalpySolver.Solve(tempEstimate);

                Temperature = new Temperature(temperature).Convert(Component.TemperatureUnits);
                Pressure = pressure.Convert(Component.PressureUnits);
            }
        }

        internal void ApplyIsentropicCompression(Pressure finalPressure)
        {
            //using (var managedEntropySolver = Expression.ObjectPoolContainer.Acquire<EntropyConservationLawTemperatureSolver>())
            //{
            //    var temperatureSolver = managedEntropySolver.Unwrap();
            //    temperatureSolver.Function.Configure(this, BaseTemperature.Value, BasePressure.Value, finalPressure);
            //    double temperature = temperatureSolver.Solve(BaseTemperature.Value);
            //    double compressionRatio = (finalPressure * BaseTemperature.Value / BasePressure.Value / temperature);
            //    this.Temperature = new Temperature(temperature).Convert(Component.TemperatureUnits);
            //    this.Pressure = new Pressure(finalPressure).Convert(Component.PressureUnits);
            //    this.CompressionRatio = new CompressionRatio(compressionRatio);
            //}

            using (var solverObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionSoneSolver>())
            using( var vectorObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionVector>())
            {
                var solver = solverObject.Unwrap();
                var vector = vectorObject.Unwrap();
                var ratios = solver.Function.Configure(this, finalPressure);
                vector.CompressionRatio = ratios.DensityRatio;
                vector.Pressure = finalPressure;
                vector.Temperature = BaseTemperature.Value * ratios.TemperatureRatio;
                object anchor = null;
                solver.Solve(vector, ref anchor);
                this.Temperature = new Temperature(vector.Temperature).Convert(Component.TemperatureUnits);
                this.Pressure = finalPressure.Convert(Component.PressureUnits);
                this.CompressionRatio = new CompressionRatio( vector.CompressionRatio);
            }

        }

        internal void ApplyIsentropicCompression(Temperature temperature)
        {
            //double ds0 = CalcEntropy(temperature) - CalcEntropy(BaseTemperature.Value);
            //double pressure = BasePressure.Value * Math.Exp(ds0 / Constants.Rgas);
            //double compressionRatio = (pressure * BaseTemperature.Value) / (BasePressure.Value * temperature);
            //this.Temperature = new Temperature(temperature).Convert(Component.TemperatureUnits);
            //this.Pressure = new Pressure(pressure).Convert(Component.PressureUnits);
            //this.CompressionRatio = new CompressionRatio(compressionRatio);

            using (var solverObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionSoneSolver>())
            using (var vectorObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionVector>())
            {
                var solver = solverObject.Unwrap();
                var vector = vectorObject.Unwrap();
                var ratios = solver.Function.Configure(this, temperature);
                vector.CompressionRatio = ratios.DensityRatio;
                vector.Pressure = BasePressure.Value * ratios.PressureRatio;
                vector.Temperature = temperature;
                object anchor = null;
                solver.Solve(vector, ref anchor);
                this.Temperature = new Temperature(vector.Temperature).Convert(Component.TemperatureUnits);
                this.Pressure = new Pressure(vector.Pressure).Convert(Component.PressureUnits);
                this.CompressionRatio = new CompressionRatio(vector.CompressionRatio);
            }
        }

        internal void ApplyIsentropicCompression(CompressionRatio compressionRatio)
        {
            //using (var managedEntropySolver = Expression.ObjectPoolContainer.Acquire<EntropyConservationLawTemperatureSolver>())
            //{
            //    var temperatureSolver = managedEntropySolver.Unwrap();
            //    temperatureSolver.Function.Configure(this, BaseTemperature.Value, BasePressure.Value, compressionRatio);
            //    double temperature = temperatureSolver.Solve(BaseTemperature.Value);
            //    double pressure = compressionRatio * BasePressure.Value * temperature / BaseTemperature.Value;
            //    this.Temperature = new Temperature(temperature).Convert(Component.TemperatureUnits);
            //    this.Pressure = new Pressure(pressure).Convert(Component.PressureUnits);
            //}

            using (var solverObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionSoneSolver>())
            using (var vectorObject = Expression.ObjectPoolContainer.Acquire<IsentropicCompressionVector>())
            {
                var solver = solverObject.Unwrap();
                var vector = vectorObject.Unwrap();
                var ratios = solver.Function.Configure(this, compressionRatio);
                vector.CompressionRatio = compressionRatio;
                vector.Pressure = BasePressure.Value * ratios.PressureRatio;
                vector.Temperature = BaseTemperature.Value * ratios.TemperatureRatio;
                object anchor = null;
                solver.Solve(vector, ref anchor);
                this.Temperature = new Temperature(vector.Temperature).Convert(Component.TemperatureUnits);
                this.Pressure = new Pressure(vector.Pressure).Convert(Component.PressureUnits);
                this.CompressionRatio = new CompressionRatio(vector.CompressionRatio);
            }
        }



        internal double CalcDensity(double temperature, double pressure)
        {
            double molarConcentration = pressure / (Constants.Rgas * temperature);
            return molarConcentration * this.MolarMass;
        }

        internal double CalcCp(double temperature)
        {
            using (var arrayObject = Expression.SpeciesArrayObjectPool.Acquire())
            {
                double[] specCp = arrayObject.Unwrap();
                Expression.CalcCpHeatCapacities(specCp, temperature);
                double cp = 0.0;
                for (int isp = 0; isp < Expression.Nsp; ++isp)
                {
                    cp += m_spX[isp] * specCp[isp];
                }
                return cp * Constants.Rgas;
            }
        }

        internal double CalcGamma(double temperature)
        {
            double cp = CalcCp(temperature);
            double gamma = cp / (cp - Constants.Rgas);
            return gamma;
        }

        internal double CalcSoundVelocity(double temperature)
        {
            double gamma = CalcGamma(temperature);
            return Math.Sqrt( gamma * Constants.Rgas * temperature / this.MolarMass );
        }



        internal void ApplyShockWaveCompression(Temperature? temperature, Pressure? pressure)
        {
            using (var poolItemSolver = Expression.ObjectPoolContainer.Acquire<ShockWaveSoneSolver>())
            using (var poolItemVector = Expression.ObjectPoolContainer.Acquire<ShockCompressionVector>())
            {
                var solver = poolItemSolver.Unwrap();
                var vector = poolItemVector.Unwrap();

                double gamma = this.CalcGamma(temperature.HasValue ? temperature.Value : this.BaseTemperature.Value);

                vector.MachNumber = this.MachNumber.HasValue ? this.MachNumber.Value : 2.0;
                vector.CompressionRatio = this.CompressionRatio.HasValue ? this.CompressionRatio.Value : 50.0;
                vector.Temperature = temperature.HasValue ? temperature.Value.ToDouble() : (5000.0 + this.BaseTemperature.Value.ToDouble());
                vector.Pressure = pressure.HasValue ? pressure.Value.ToDouble() : (10.0 * this.BasePressure.Value.ToDouble());

                if (this.MachNumber.HasValue)
                {
                    if (MachNumber.Value < 1.0) { throw new InvalidOperationException(); }
                    if (MachNumber.Value < 1.00001)
                    {
                        CompressionRatio = new CompressionRatio(1.0);
                        Temperature = BaseTemperature.Value;
                        Pressure = BasePressure.Value;
                        return;
                    }
                    var ratios = PolytropicGasRatios.CaclShockWaveRatios(this.MachNumber.Value, gamma);
                    vector.MachNumber = this.MachNumber.Value;
                    vector.Pressure = this.BasePressure.Value * ratios.PressureRatio;
                    vector.Temperature = this.BaseTemperature.Value * ratios.TemperatureRatio;
                    vector.CompressionRatio = ratios.DensityRatio;
                }
                else if (this.CompressionRatio.HasValue)
                {
                    if (CompressionRatio.Value < 1.00001)
                    {
                        MachNumber = new MachNumber(1.0);
                        Temperature = BaseTemperature.Value;
                        Pressure = BasePressure.Value;
                        return;
                    }
                    MachNumber m = PolytropicGasRatios.CalcMachNumber(this.CompressionRatio.Value, gamma);
                    if (double.IsNaN(m)) { m = new MachNumber(10.0); }
                    var ratios = PolytropicGasRatios.CaclShockWaveRatios(m, gamma);
                    vector.MachNumber = m;
                    vector.Pressure = this.BasePressure.Value * ratios.PressureRatio;
                    vector.Temperature = this.BaseTemperature.Value * ratios.TemperatureRatio;
                    vector.CompressionRatio = this.CompressionRatio.Value;
                }
                else if (temperature.HasValue)
                {
                    DoubleRelativeTolerance rt = new DoubleRelativeTolerance(1.0e-5);
                    if (true == rt.AreSame(temperature.Value, BaseTemperature.Value))
                    {
                        MachNumber = new MachNumber(1.0);
                        this.CompressionRatio = new Nature.CompressionRatio(1.0);
                        Temperature = temperature.Value;
                        Pressure = BasePressure.Value;
                        return;
                    }
                    var m = PolytropicGasRatios.CalcMachNumber(BaseTemperature.Value, temperature.Value, gamma);
                    var ratios = PolytropicGasRatios.CaclShockWaveRatios(m, gamma);
                    vector.Temperature = temperature.Value;
                    vector.MachNumber = m;
                    vector.CompressionRatio = ratios.DensityRatio;
                    vector.Pressure = BasePressure.Value * ratios.PressureRatio;
                }
                else if (pressure.HasValue)
                {
                    DoubleRelativeTolerance rt = new DoubleRelativeTolerance(1.0e-5);
                    if (true == rt.AreSame(pressure.Value, BasePressure.Value))
                    {
                        MachNumber = new MachNumber(1.0);
                        this.CompressionRatio = new Nature.CompressionRatio(1.0);
                        Temperature = BaseTemperature.Value;
                        Pressure = pressure.Value;
                        return;
                    }
                    var m = PolytropicGasRatios.CalcMachNumber(BasePressure.Value, pressure.Value, gamma);
                    var ratios = PolytropicGasRatios.CaclShockWaveRatios(m, gamma);
                    vector.Pressure = pressure.Value;
                    vector.MachNumber = m;
                    vector.CompressionRatio = ratios.DensityRatio;
                    vector.Temperature = BaseTemperature.Value * ratios.TemperatureRatio;
                }

                solver.Function.Configure(this, temperature, pressure);
                object anchor = null;
                solver.Solve(vector, ref anchor);
                this.MachNumber = new MachNumber(vector.MachNumber);
                this.CompressionRatio = new CompressionRatio(vector.CompressionRatio);
                this.Temperature = new Temperature(vector.Temperature).Convert(Component.TemperatureUnits);
                this.Pressure = new Pressure(vector.Pressure).Convert(Component.PressureUnits);
            }
        }

    }
}
    